function [res, fout] = fit_func(lam,x,y,model,alpha,Vr)
% fit_func switches between fitting models

switch (model)
    case 1
        %if model==1
        % SINGLE exponentials
        %	assumes a function of the form
        %	  y =  LAM1(1)*exp(-lam(2)*x)
        f=lam(1)*exp(-lam(2)*x);
        res = y-f;
    case 2
        % if model==2 double exp
        %	with 2 linear parameters 2 nonlinear parameters .
        f=lam(1).*exp(-x/lam(2))+lam(3).*exp(-x/lam(4));
        res = y-f;
        
    case 3
        % model for fitting hyperpolarizing voltage steps at misc times
        % 1/7/99 P. Manis Model has 2 exp decay plus offset, and value =
        % offset (lam1) when x < alpha (t0) else value  = lam1 - 2 exp
        % decays (amplitudes adjustable)
        %
        x1 = x(find(x < alpha)); %#ok<FNDSB>
        z1 = lam(1)+0*x1; % make same length as z1
        x2 = x(find(x >= alpha)); %#ok<FNDSB> % the rest of passed points
        z2 = lam(1) + lam(2).*(1-exp(-(x2-alpha)/lam(3))) + ...
            lam(4).*(1-exp(-(x2-alpha)/lam(5)))+lam(6);
        f = cat(1,z1,z2);   % combine piecwise.
        res = y-f;
        
    case 4
        % model for fitting hyperpolarizing voltage steps with single
        % exponential...
        x1 = x(find(x < alpha)); %#ok<FNDSB>
        z1 = lam(1)+0*x1; % make same length as z1
        x2 = x(find(x >= alpha)); %#ok<FNDSB> % should be the rest of passed points
        z2 = lam(1) + lam(2).*(1-exp(-(x2-alpha)/lam(3)))+lam(4);
        f = cat(1,z1,z2);   % combine piecwise.
        res = y-f;



    case 9
        %if model==9
        % multiple exponentials
        %	assumes a function of the form
        %	  y =  c(1)*exp(-lam(1)*t) + ... + c(n)*exp(-lam(n)*t)
        %	with n linear parameters c and n nonlinear parameters lam.
%         n=length(lam);
        A = zeros(length(t),length(lam));
        for j = 1:size(lam)
            A(:,j) = exp(-lam(j)*x);
        end
        c = A\y;
        z = A*c;
        f = z;
        res = y-f;




        % Boltzmans in the next sectoin alpha=temperature factor (0.04 for
        % 22 C)

    case 10
        %if model==10
        % normalized  Boltzman lam(1)=z lam(2)=V0.5
        f=1./(1+exp(-lam(1)*alpha*(x-lam(2))));
            res = y-f;
    case 11
        % if model==11 normalized  Boltzman corrected for reversal
        % potential lam(1)=z lam(2)=V0.5
        f=1./(1+exp(-lam(1)*alpha*(x-lam(2)))).*(x-Vr);
         res = y-f;
   case 12
        %if model==12
        % normalized  Boltzman lam(1)=z lam(2)=V0.5
        %lam(3) =offset
        f=(1./(1+exp(-lam(1)*alpha*(x-lam(2))))+lam(3));
        res = y-f;

    case 20
        % if model==20 non normalized  Boltzman
        %lam(1)=Imax
        % lam(2)=z lam(2)=V0.5
        f=lam(1)./(1+exp(-lam(2)*alpha*(x-lam(3))));
        res = y-f;

    case 21
        %if model==21 lam(1)=Imax
        % non normalized  Boltzman in current space, with reversal
        % potential lam(2)=z lam(3)=V0.5
        f=(lam(1).*(x-Vr))./(1+exp(-lam(2)*alpha*(x-lam(3))));
        res = y-f;

    case 22
        %if model==22
        % non normalized  Boltzman with offset
        %lam(1)=Imax
        % lam(2)=z lam(4)=V0.5
        %lam(4)= offset
        f=(lam(4)+lam(1)./(1+exp(-lam(2)*alpha*(x-lam(3)))));
        res = y-f;

    case 23
        %if model==23 lam(1)=Imax
        % non normalized  Boltzman corrected for reversal potential
        % lam(2)=z lam(3)=V0.5 lam(4)=offset
        f=(lam(1)./(1+exp(-lam(2)*alpha*(x-lam(3)))).*(x-Vr)+lam(4));
        res = y-f;

    case 24
        % if model==24
        %lam(1)=Imax
        % non normalized  Boltzman raised to power 4 corrected for reversal
        % potential lam(2)=z lam(3)=V0.5 lam(4)=offset
        f=(lam(1)./((1+exp(-lam(2)*alpha*(x-lam(3)))).^4).*(x-Vr)+lam(4));
        res = y-f;

    case 25
        %if model==25
        % non normalized  Double Boltzman with offset
        %lam(1)=Imax
        % lam(2)=z lam(4)=V0.5
        %lam(4)= offset
        f=(lam(7)+(lam(1)./(1+exp(-lam(2)*alpha*(x-lam(3))))) + ...
            (lam(4)./(1+exp(-lam(5)*alpha*(x-lam(6))))));
        res = y-f;

    case 30
        % if model==30 rising single exponential
        f=(lam(1)*(1-exp(-lam(2)*x)));
         res = y-f;

    case 40
        % if model==40
        %   Boltzman
        % lam(1)=slope lam(2)=V0.5
        f=1./(1+exp(-(x-lam(2))/lam(1)));
        res = y-f;

    case 41
        % if model==41 normalized  Boltzman corrected for reversal
        % potential lam(1)=slope lam(2)=V0.5
        f=1./(1+exp(-(x-lam(2))/lam(1))).*(x-Vr);
            res = y-f;

    case 42
        % if model==42
        %   Boltzman
        % lam(1)=slope lam(2)=V0.5
        %lam(3)=offset
        f=(1./(1+exp(-(x-lam(2))/lam(1)))+lam(3));
            res = y-f;

    case 44
        % if model==44 normalized  Boltzman raied to the 4th power
        % corrected for reversal potential lam(1)=slope lam(2)=V0.5
        f=(1./((1+exp(-(x-lam(2))/lam(1)))).^4).*(x-Vr);
            res = y-f;

    case 45
        % function for HWF data anlaysis 1/6/98 P. Manis Boltzman expressed
        % in reversed coordinates for fitting to current-clamp IV includes
        % a linear term. lam(1) is gmax, lam(3) is vhalf, lam(2) is slope
        % lam(4) is gleak, lam(5) is vleak.
        f=(((lam(1).*(x-Vr))./(1+exp(-(x-lam(3))*lam(2)*alpha)))+ ...
            (lam(4).*(x-lam(5))));
        res = y-f;

    case 50
        % if model==50 e single exponential decay
        f=(lam(1)+lam(2).*exp(-x/lam(3)));
        res = y-f;

    case 51
        % if model==51 double exponential decay
        f=lam(1)+lam(2).*exp(-x/lam(3))+lam(5).*exp(-x/lam(4));
        res = y-f;

    case 52
        % if model==52 exponential acivation and single exponential decay;
        % activation has power
        f = lam(4) + lam(1) .* ((1-exp(-x/lam(2))).^alpha) .* ...
            exp(-x/lam(3));
        res = y-f;

    case 53
        % if model==53 single exponential activation and double exponential
        % decay activation has power alpha (fixed).
        f=lam(6)+lam(1).*((1-exp(-x/lam(2))).^alpha).*((1-lam(5))* ...
            exp(-x/lam(3))+(lam(5)*exp(-x/lam(4))));
        res = y-f;

    case 54 % (this is the same function as 53...)
        % if model==54 PBM EPEE Model used to fit A currents
        f=(lam(1)+lam(2).*((1-exp(-x/lam(3))).^alpha).*(lam(4).* ...
            (exp(-x/lam(5)))+(1-lam(4)).*exp(-x/lam(6))));
        res = y-f;


    case 55
        % if model==55
        %  double exponential decay
        f=lam(1).*(1-(lam(2).*exp(-x/(lam(3)/1000)) + ...
            lam(4).*exp(-x/(lam(5)/1000))));
                res = y-f;
% end
    case 56
        % if model==56
        %  double exponential decay + const
        f=(lam(1)+(lam(2).*(lam(3).*exp(-x/lam(4)) + ...
            (1-lam(3)).*exp(-x/lam(5)) ) ));
                res = y-f;
% end

        % PBM Modified EPEE Model used to fit A currents consists of m^k*h
        % + n where n is first order activation of slower ss component...
        %
    case 57
        % if model==57
        f=(lam(1)+lam(2).*((1-exp(-x/lam(3))).^alpha).*(exp(-x/lam(4)))+lam(5).*(1-exp(-x/lam(6))));
                res = y-f;
% end
        % PBM Modified EPEE Model used to fit A currents consists of
        % m^k*(1-p*h) where h is first order inactivation and has a
        % non-inactivating fraction determined by p
        %
    case 58
        % if model==58
        f=(lam(1)+(lam(2).*((1-exp(-x/lam(3))).^alpha).*(lam(5).*(exp(-x/lam(4)))+(1-lam(5)))));
                res = y-f;
% end
    case 61
        % if model==61 double exponential decay, with time constants fixed
        f=(lam(1)+lam(2).*exp(-x/alpha)+lam(3).*exp(-x/Vr));
                res = y-f;
% end
    case 62
        % if model==62 single exponential decay, with time constants fixed
        % fit 3 traces simultaneously
        f1=y(1,:)-(lam(1)+lam(2).*exp(-x(1,:)/lam(3)));
        f2=y(2,:)-(lam(4)+lam(5).*exp(-x(2,:)/lam(3)));
        f3=y(3,:)-(lam(6)+lam(7).*exp(-x(3,:)/lam(3)));
        dd=1.5;
        if lam(1) < 0 || lam(2) < 0
            f1=f1.*(dd);
        end
        if lam(4) < 0 || lam(5) < 0
            f2=f2.*(dd);
        end
        if lam(6) < 0 || lam(7) < 00
            f3=f3.*(dd);
        end
  %      f=(f1+f2+f3);
        f=(f1.^2+f2.^2+f3.^2).^(0.5);
        res = f;
        % end
    case 63
        % if model==63 double exponential decay, with time constants fixed
        % fit 3 traces simultaneously
        f1=y(1,:)-(lam(1)+lam(2).*exp(-x(1,:)/lam(3))+lam(4).*exp(-x(1,:)/lam(5)));
        f2=y(2,:)-(lam(6)+lam(7).*exp(-x(2,:)/lam(3))+lam(8).*exp(-x(2,:)/lam(5)));
        f3=y(3,:)-(lam(9)+lam(10).*exp(-x(3,:)/lam(3))+lam(11).*exp(-x(3,:)/lam(5)));
        dd=1.5;
        if lam(1) < 0 || lam(2) < 0 || lam(4) < 0
            f1=f1.*(dd);
        end
        if lam(6) < 0 || lam(7) < 0 || lam(8) < 0
            f2=f2.*(dd);
        end
        if lam(9) < 0 || lam(10) < 0 || lam(11) < 0
            f3=f3.*(dd);
        end
        % f=(f1+f2+f3);
        f=(f1.^2+f2.^2+f3.^2).^(0.5);
         res = f;


    case 72
        % if model==72 single exponential decay, with time constants fixed
        % fit 6 traces simultaneously
        f1=y(1,:)-(lam(1)+lam(2).*exp(-x(1,:)/lam(3)));
        f2=y(2,:)-(lam(4)+lam(5).*exp(-x(2,:)/lam(3)));
        f3=y(3,:)-(lam(6)+lam(7).*exp(-x(3,:)/lam(3)));
        f4=y(4,:)-(lam(8)+lam(9).*exp(-x(4,:)/lam(3)));
        f5=y(5,:)-(lam(10)+lam(11).*exp(-x(5,:)/lam(3)));
        f6=y(6,:)-(lam(12)+lam(13).*exp(-x(6,:)/lam(3)));
        dd=1.5;
        if lam(1) < 0 || lam(2) < 0
            f1=f1.*(dd);
        end
        if lam(4) < 0 || lam(5) < 0
            f2=f2.*(dd);
        end
        if lam(6) < 0 || lam(7) < 00
            f3=f3.*(dd);
        end
        if lam(8) < 0 || lam(9) < 00
            f4=f4.*(dd);
        end
        if lam(10) < 0 || lam(11) < 00
            f5=f5.*(dd);
        end
        if lam(12) < 0 || lam(13) < 00
            f6=f6.*(dd);
        end
        f=(f1.^2+f2.^2+f3.^2+f4.^2+f5.^2+f6.^2).^(0.5);
                res = f;
% end
    case 73
        % if model==73 double exponential decay, with time constants fixed
        % fit 6 traces simultaneously
        f1=y(1,:)-(lam(1)+lam(2).*exp(-x(1,:)/lam(3))+lam(4).*exp(-x(1,:)/lam(5)));
        f2=y(2,:)-(lam(6)+lam(7).*exp(-x(2,:)/lam(3))+lam(8).*exp(-x(2,:)/lam(5)));
        f3=y(3,:)-(lam(9)+lam(10).*exp(-x(3,:)/lam(3))+lam(11).*exp(-x(3,:)/lam(5)));
        f4=y(4,:)-(lam(12)+lam(13).*exp(-x(4,:)/lam(3))+lam(14).*exp(-x(4,:)/lam(5)));
        f5=y(5,:)-(lam(15)+lam(16).*exp(-x(5,:)/lam(3))+lam(17).*exp(-x(5,:)/lam(5)));
        f6=y(6,:)-(lam(18)+lam(19).*exp(-x(6,:)/lam(3))+lam(20).*exp(-x(6,:)/lam(5)));
        dd=2;
        if lam(1) < 0 || lam(2) < 0 || lam(4) < 0
            f1=f1.*(dd);
        end
        if lam(6) < 0 || lam(7) < 0 || lam(8) < 0
            f2=f2.*(dd);
        end
        if lam(9) < 0 || lam(10) < 0 || lam(11) < 0
            f3=f3.*(dd);
        end
        if lam(12) < 0 || lam(13) < 0 || lam(14) < 0
            f4=f4.*(dd);
        end
        if lam(15) < 0 || lam(16) < 0 || lam(17) < 0
            f5=f5.*(dd);
        end
        if lam(18) < 0 || lam(19) < 0 || lam(20) < 0
            f6=f6.*(dd);
        end
        f=(f1.^2+f2.^2+f3.^2+f4.^2+f5.^2+f6.^2).^(0.5);
                res =f;
% end

    case 74
        % if model==74 double exponential decay, with time constants fixed
        % fit N traces simultaneously lam(1) is tau1 (SHARED) lam(2) is
        % tau2 (SHARED) lam(3, 4, 5) are offset, a1 and a2 for trace 1
        % lam(6, 7, 8) are offsest, a1 and a2 for trace 2 etc.
        f = 0;
        for i = 1:size(y,1)
            j = i*3;
            fn = y(i,:) - (lam(j)+lam(j+1).*exp(-x(i,:)*lam(1))+lam(j+2).*exp(-x(i,:))*lam(2));
            f = f + (fn.^2);
        end;
        f = f.^0.5;
        res = f;

    case 75
        % if model == 75, Woodhull fit from Washburn et al, J. Neurosci,
        % 1998. input is IV for current parameters are: lam(1) is g scaling
        % factor lam(2) is blocker concentration lam(3) is Kd(0) lam(4) is
        % z*(1-delta) alpha is RT/F (must be fixed when called!) Vr is
        % reversal potential... G0 is set to be the function given in that
        % paper...
        g0 = 0.3 + 0.6*exp((x-50)/40);
        f = lam(1).*(x-Vr).*g0./(1 + (lam(2)./(lam(3)*exp(-lam(4)*x/alpha))));
        res = y-f;

    case 76
        % THIS MODEL is IN PLANING (P. Manis, 2006). if model==76 double
        % exponential decay, with time constants fixed amplitudes must vary
        % according to a Boltzmann function (e.g., parameterized
        %
        % fit N traces simultaneously lam(1) is tau1 (SHARED) because at
        % fixed voltage lam(2) is tau2 (SHARED) because at fixed voltage
        % lam(3) is A1 (shared boltzman max1) lam(4) is I1 (shared boltzman
        % offset1) lam(5) is V1 (shared boltzman vhalf1) lam(6) is Z1
        % (shared boltzman slope1) lam(7) is A1 (shared boltzman max2)
        % lam(8) is I1 (shared boltzman offset2) lam(9) is V1 (shared
        % boltzman vhalf2) lam(10) is Z1 (shared boltzman slope2) lam(5, 6,
        % 7, ... ) are the offsets for each trace 1 Requires passing z (V
        % for the boltzman to work with on each trace). etc.
        u = size(y);
        f = 0;
        for i = 1:u(1)
            j = i*3;
            fn = y(i,:) - (lam(j)+((lam(j+1).*exp(-x(i,:)/lam(1)))+((lam(j+2)).*exp(-x(i,:))/lam(2))));
            %      f = f + (fn.^2); % compute sum of squares in error
            f = f + (fn.^2);
        end;
        f = f.^0.5;
        res = f;
    otherwise
end
if(nargout == 2)
    fout = f;
end;
return;

